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Abstract.  Hydrogen  plumes  gained  much  interest  from  the  fact  that  H2  is  an  important  reaction  product 
of  a  hydrazine  thruster.  H2  shows  strong  and  temperature  dependant  rotational  relaxation  which  quickly 
leads  to  nonequilibrium  between  translational  and  rotational  temperature  in  expanding  flows.  In  a  hybrid 
calculation  of  the  whole  flowfield,  both  methods  (Navier-Stokes  and  DSMC)  must  model  rotational  relaxation 
consistently.  Jeans’  equation  was  used  for  Navier-Stokes,  and  for  DSMC  the  Borgnakke-Larsen  model  with 
constant  and  with  energy-dependent  relaxation  probability  was  used.  The  comparison  with  free  molecule 
probe  measurements  of  hydrogen  plumes  in  the  DLR  High  Vacuum  Facility  STG  provides  insight  into  the 
accuracy  of  the  numerical  models. 


EXPERIMENTAL  SET  UP 

Hydrogen  is  an  important  species  in  the  exhaust  plume  of  a  hydrazine  thruster.  To  prepare  tests  with  real 
hydrazine  thrusters,  simulated  plumes  of  pure  hydrogen  have  been  studied  in  the  DLR  High  Vacuum  Facility 
STG  [1]  (see  fig.  1).  The  very  high  vacuum  (<  10-5  mbar)  provided  by  the  STG  facility  during  thruster 
operation  enables  measurements  with  a  free  molecule  pressure  probe  with  negligible  background  disturbances 
throughout  the  whole  flowfield.  Because  of  the  low  boiling  point  of  H2,  the  cryo-pump  which  surrounds  the 
expansion  room  has  to  be  operated  with  liquid  Helium. 


FIGURE  1.  Cross  section  of  the  STG  vacuum  chamber  with  cold  walls,  thruster,  probe  and  linear  positioners. 


The  conical  0.5  N  nozzle  and  the  free  molecule  pressure  probe  can  be  rotated  and  the  probe  can  also  be  moved 
in  radial  direction  from  the  nozzle.  Angular  profiles  of  number  flux  were  obtained  by  rotating  the  nozzle  and 
radial  profiles  at  different  angles  by  moving  the  probe.  The  slit  of  the  probe  could  have  different  orientations 
relative  to  the  nozzle,  but  here  the  probe  was  always  directly  facing  the  nozzle.  Stagnation  pressure  was  held  at 
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Po  =  0.95  bar  and  temperature  To  =  300JRC.  These  conditions  yield  a  nozzle  exit  Reynolds  number  of  Re  =  800. 
The  nozzle  wall  had  a  temperature  of  300if. 


CONTINUUM  METHOD 

The  hydrogen  flow  inside  the  nozzle  was  calculated  with  the  DLR  Navier-Stokes  solver  CEVCATS,  a  block- 
structured  3D  upwind  finite-volume  code  developed  at  DLR  [2,3].  CEVCATS  uses  a  hybrid  AUSM  -  van 
Leer  method  for  flux  calculation  and  multigrid  methods  with  residual  averaging  and  local  time-stepping  to 
accelerate  convergence  to  steady  state.  The  modelling  of  the  rotational  energy  relaxation  was  implemented 
following  Jeans’  equation  [4], 


deT 

dt 


tt  —  ZTrc, 


(1) 


where  e*(t)  is  the  instantaneous  equilibrium  rotational  energy  which  is  e*(t)  =  er(oo)  =  kTt  in  the  isothermal 
limit,  rc  is  the  instantaneous  collision  time.  The  rotational  collision  number  Zr  is  either  constant  or  a  function 
of  the  translational  temperature  T.  The  temperature  and  viscosity  coefficient  w  dependant  Zr  is  modelled  as 
in  the  paper  of  Boyd  [5]  where  good  agreement  with  experimental  data  cited  by  Lambert  [6]  was  shown: 


ZT  =  10480/T“., 


(2) 


INTERFACING  OF  THE  TWO  METHODS 

Parts  of  the  boundary  layer  inside  the  nozzle  and  the  flowfield  outside  the  nozzle  where  the  assumptions 
for  the  Navier-Stokes  equations  cease  to  be  valid  have  to  be  calculated  with  DSMC.  The  exact  location  of 
this  interface  is  determined  by  a  parameter  characterizing  the  deviation  of  the  distribution  function  /  from  a 
Maxwell  distribution  [7,8].  A  parameter  which  characterizes  the  breakdown  of  continuum  is  the  maximum  of 
the  viscous-stress  tensor  ntj  and  the  heat  flux  vector  qp. 

B  =  max^Tij],  |<7i|).  (3) 

The  interface  location  was  choosen  such  that  this  parameter  didn’t  exceed  0.1.  Figure  2  shows  B  along  the 
interface  which  was  choosen  as  a  straight  line  from  the  nozzle  axis  to  the  nozzle  wall. 


y  [mm] 

FIGURE  2.  Continuum  breakdown  parameters  along  the  method  interface 
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B  stays  very  low  in  the  core  flow  but  increases  rapidly  in  the  boundary  layer  because  of  the  high  viscous 
stress.  Towards  the  nozzle  wall,  the  temperature  gradient  and  therefore  the  heat  flux  becomes  dominant.  As 


a  reference,  the  local  Knudsen  number  Kni„cai  =  ~  ^  [7]  is  also  plotted.  Kniocai  shows  the  same  increase  in 
the  boundary  layer  but  doesn’t  predict  the  nonequilibrium  conditions  near  the  nozzle  wall  due  to  the  strong 
temperature  gradient. 


DSMC  CALCULATION 

The  DSMC  part  of  the  flowfield  was  computed  with  the  MONACO  [9]  code.  The  Natural-Sample-Size 
method  by  Baganoff  and  McDonald  [10]  is  used  as  collision  sampling  technique.  Collisions  are  simulated  by 
the  Variable  Hard  Sphere  (VHS)  and  the  Larsen-Borgnakke  model  [11].  A  simple  diffuse  reflection  model  with 
full  thermal  accommodation  handles  gas-surface  interaction.  The  code  runs  on  different  parallel  architectures 
like  Hitachi  SR8000  and  Sun  Sparc  MP. 

The  unstructured  grid  [12]  was  automatically  adapted  to  the  local  mean  free  path,  resulting  in  a  grid  size  of 
92000  cells.  The  flowfield  was  simulated  with  7.7  million  particles  and  one  calculation  took  90  h  of  computing 
time  on  one  node  of  a  Hitachi  SR8000. 

The  inflow  conditions  for  the  DSMC  part  are  given  by  the  flow  conditions  of  the  Navier-Stokes  solution  at  the 
interface.  Here  the  flow  is  clearly  not  equilibrium  any  more,  so  the  inflowing  particles  have  to  be  generated  from 
a  Chapman-Enskog  distribution.  An  acceptance-rejection  random  velocity  generator  for  the  Chapman-Enskog 
distribution  as  proposed  by  Garcia  and  Alder  [13]  is  used. 

Fig.  3  shows  the  number  density  as  given  by  the  Navier-Stokes  solution  and  the  number  density  at  the 
interface,  sampled  from  the  first  row  of  DSMC  cells.  The  inflowing  particles  were  generated  from  a  Maxwell 
distribution  and  from  the  Chapman-Enskog  distribution.  While  the  densities  agree  very  well  in  the  core  flow 
where  the  condition  is  near  equilibrium,  the  particles  generated  from  a  Maxwell  distribution  lead  to  a  much 
lower  density  near  the  nozzle  wall  than  those  generated  from  a  Chapman-Enskog  distribution  which  reproduce 
the  density  of  the  Navier-Stokes  solution.  This  clearly  shows  that  for  this  case,  inflowing  particles  should  be 
generated  from  a  Chapman-Enskog  distribution  to  model  the  inflow  conditions  correctly. 
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FIGURE  3.  Number  densities  at  the  interface  as  given  by  the  Navier-Stokes  solution  and  sampled  from  the  first  row 
of  DSMC  cells. 


It  can  be  shown  as  in  a  paper  by  Boyd  et  al.  [5]  that  an  energy-dependant  form  of  the  probability  of  rotational 
energy  exchange  reduces  to  (2)  when  integrated  over  the  equilibrium  distribution  function, 


i  r(£  +  5/2  —  u)  (EAU 
10480  r(C  +  5/2)  \k)  ’ 


(4) 


821 


where  Ec  is  the  total  collision  energy,  k  is  the  Boltzmann  constant  and  (  is  the  number  of  degrees  of  freedom 
participating  in  the  energy  exchange  process.  Due  to  the  difference  between  the  probability  of  rotational  energy 
exchange  in  DSMC  and  the  macroscopic  relaxation  rate,  some  additional  adjustments  have  to  be  made  [5,14]. 


RESULTS 


Figure  4a, b  shows  Navier-Stokes  and  DSMC  calculations  of  translational  and  rotational  temperature  with 
temperature-dependant  relaxation  probablility  in  the  vicinity  of  the  nozzle. 


FIGURE  4.  Navier-Stokes/DSMC  calculation  of  a  0.5  N  conical  nozzle  flow  with  Re®  =  800,  testgas  H2.  Contours 
of  translational  (a)  and  rotational  temperature  with  temperature  dependant  (b)  and  constant  (c)  rotational  collision 
number  are  shown. 
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Starting  near  the  nozzle  throat,  the  rotational  temperature  quickly  departs  from  the  translational  tempera¬ 
ture  as  expected  for  hydrogen.  Figure  4c  shows  the  rotational  temperature  calculated  with  constant  ZT  =  500. 
Compared  to  4b,  the  relaxation  is  faster  near  the  nozzle  throat  and  slower  after  the  nozzle  exit.  Although  the 
number  of  simulated  molecules  was  too  low  in  the  backflow  region  to  produce  reliable  rotational  temperatures, 
the  number  density  in  this  region  indicates  that  for  constant  ZT,  more  molecules  are  scattered  at  higher  angles 
into  the  backflow  due  to  the  delayed  transfer  of  rotational  energy  into  the  translational  modes. 

The  comparison  between  measured  and  calculated  number  flux  on  radial  and  angular  profiles  is  shown 
in  figure  5.  The  total  number  flux,  obtained  by  integrating  the  angular  profiles  over  a  sphere,  shows  good 
agreement  between  measurements  and  calculation.  However,  the  plume  shape  as  shown  in  the  angular  profile 
is  predicted  less  slender  by  the  calculation  than  the  measurements  indicate.  This  leads  to  a  lower  calculated 
number  flux  near  the  axis  and  somewhat  higher  flux  for  higher  angles.  The  asymmetry  in  the  measured 
plume  comes  from  the  influence  of  the  high  density  measurements  near  the  axis  on  the  later  measurements  at 
larger  negative  angles  due  to  incomplete  outgassing.  The  radial  profiles  show  the  same  overprediction  of  the 
calculation  at  higher  angles.  For  the  90°  case,  the  measurements  show  no  decrease  in  number  flux  for  higher 
radial  distances  which  again  indicates  incomplete  outgassing  of  the  probe. 


Radial  distance  [mm]  Angie  theta  [degl 

FIGURE  5.  Patterson  probe  measuremens  and  calculations  of  i?2  number  flux  on  radial  and  angular  profiles. 


Variation  of  the  model  parameters  for  rotational  relaxation  showed  no  significant  influence,  when  compared 
with  the  number  flux  measurements  in  the  available  singular  range.  To  decide  on  these  details,  direct  mea¬ 
surements  of  the  rotational  temperature  would  give  more  insight.  One  reason  for  the  differences  in  the  plume 
shape  may  be  the  nozzle  wall  temperature  which  is  only  known  at  one  point  in  the  middle  of  the  nozzle  wall 
from  the  experiments.  In  the  calculations,  the  nozzle  wall  was  assumed  to  have  a  constant  temperature  which 
may  not  reflect  the  actual  temperature  distribution. 


CONCLUSIONS 

The  hydrogen  flowfield  of  a  0.5  N  thruster  nozzle  was  calculated  with  a  decoupled  hybrid  Navier-Stokes  - 
DSMC  method.  Consistent  models  for  the  relaxation  of  rotational  temperature  were  implemented  in  both  flow 
solvers  with  temperature-dependant  and  constant  rotational  relaxation  number  Zr.  To  decide  on  the  model 
and  its  parameters,  comparison  with  measurements  of  rotational  temperature  are  needed  because  the  effect  on 
the  number  flux,  for  which  experimental  data  are  already  available,  is  too  small.  Experimental  results  of  Troi, 
density  and  velocity  measurements  with  REMPI  (Resonance  Enhanced  Multi  Photon  Ionisazion)  for  hydrogen 
will  be  available  in  the  near  future  so  the  calculations  can  be  compared  with  an  additional  independant  set  of 


823 


measurements.  It  was  shown  that  it  is  important  to  generate  particles  which  enter  the  DSMC  region  from  a 
Chapman-Enskog  distribution  because  of  the  starting  deviation  from  equilibrium  in  the  Navier-Stokes  solution. 
The  calculations  reproduced  the  total  number  flux  found  by  integrating  the  patterson  probe  measurements. 
The  plume  shape  was  in  quite  good  agreement  with  the  measurements  but  the  calculations  predicted  lower 
number  flux  on  the  axis  and  slightly  increased  number  flux  in  the  off-axis  regions. 
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